AIAA 2003-0072 


Implementation of Preconditioned Dual-Time Procedures in 

OVERFLOW 

Shishir A. Pandya * 

NASA Ames Research Center 
Moffett Field, CA 94035 

Sankaran Venkateswaran f 

University of Tennessee 
Tullahoma, TN 37388 

Thomas H. Pulliam * 

NASA Ames Research Center 
Moffett Field, CA 94035 


Abstract 

Preconditioning methods have become the method 
of choice for the solution of flowfields involving the 
simultaneous presence of low Mach and transonic re- 
gions. It is well known that these methods are impor- 
tant for insuring accurate numerical discretization as 
well as convergence efficiency over various operating 
conditions such as low Mach number, low Reynolds 
number and high Strouhal numbers. For unsteady 
problems, the preconditioning is introduced within 
a dual-time framework wherein the physical time- 
derivatives are used to march the unsteady equations 
and the preconditioned time-derivatives are used for 
purposes of numerical discretization and iterative so- 
lution. In this paper, we describe the implementation 
of the preconditioned dual-time methodology in the 
OVERFLOW code. To demonstrate the performance 
of the method, we employ both simple and practical 
unsteady flowfields, including vortex propagation in 
a low Mach number flow, flowfield of an impulsively 
started plate (Stokes’ first problem) and a cylindri- 
cal jet in a low Mach number crossflow with ground 
effect. All the results demonstrate that the precondi- 
tioning algorithm is responsible for improvements to 
both numerical accuracy and convergence efficiency 
and, thereby, enables low Mach number unsteady com- 
putations to be performed at a fraction of the cost of 
traditional time-marching methods. 
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Introduction and Background 

Accurate and efficient numerical simulation of un- 
steady fluid flows is one of the remaining challenges for 
the resolution of engineering problems. For high speed 
flows j time-marching methods developed by the aero- 
dynamics’ CFD community have been very successful 
for both Euler and Navier-Stokes problems [1,2]. Like- 
wise, for incompressible flows, unsteady algorithms 
based on projection methods [3, 4] or artificial com- 
pressibility [5, 6, 7, 8] approaches have been widely 
applied. However, many engineering problems involve 
the co-existence of both compressible and incompress- 
ible flow regimes in the same flow field. For such cases, 
generalized preconditioned dual-time approaches have 
been developed which allow the same computer code 
to be applied in all speed regimes [9, 10]. 

An example of a flow field that simultaneously in- 
volves both low speed and high speed regions is the 
Harrier aircraft in near-hover (landing approach) con- 
dition [11]. In this situation, the aircraft’s forward 
velocity is approximately 0.04 Mach. However, at the 
same time, the aircraft is kept in the air by four high 
speed jet exhausts directed toward the ground. The 
interaction of the jets with the cross-wind and the 
ground effect renders this flow field extremely time- 
dependent. A purely incompressible method would be 
appropriate for the low-speed free-stream region, but 
not for the high speed jet plumes. On the other hand, 
a transonic flow method would be appropriate for the 
jets, but is inefficient for the free-stream. Thus, nu- 
merical procedures for the solution of such problems 
must be capable of simultaneously handling both tran- 
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sonic and low speed flow regimes. 

The difficulty that transonic algorithms face at low 
speeds arises due to the poor conditioning of the eigen- 
values of the tiine-marching system. This lack of con- 
ditioning leads to poor solution accuracy as well as 
poor convergence properties. Over the last decade, 
preconditioning methods have become increasingly 
popular as a means of assuring that the time-marching 
algorithm is well-conditioned in terms of both accu- 
racy and efficiency at all speeds [12, 13]. For steady 
state problems, these methods involve the introduction 
of pseudo-time derivatives in lieu of the physical time 
derivative terms in the time-marching system. For un- 
steady flows, the formulation is typically cast within 
a dual-time-stepping strategy, wherein the physical- 
time derivatives are employed to follow the physical 
transients and the pseudo-time derivatives serve as an 
iterative device. Proper definition of the pseudo-time 
derivatives serves to optimize the numerical solution 
procedure, making the performance of the algorithm 
commensurate with traditional transonic methods at 
high speeds and with incompressible methods at low 
speeds. 

In this paper, the implementation of the precon- 
ditioned dual-time algorithm in the OVERFLOW 
code [1] is discussed. The OVERFLOW code is a com- 
pressible Navier-Stokes code that uses the Chimera 
overset grid approach [14] for simulating complex-body 
configurations such as the afore- mentioned Harrier air- 
craft. The main solution algorithm in the code is the 
diagonalized approximate-factorization procedure [15]. 
The implementation of the preconditioned dual-time 
scheme in the diagonalized approximate factorization 
framework has been described by Buelow et al. [16] and 
we follow the same approach in this paper. Specifically, 
this involves a modification of the traditional diago- 
nalization procedure to include both the physical and 
preconditioned time-derivatives, thereby avoiding the 
introduction of block matrices in the implicit operator. 

The paper is organized as follows. We start by pre- 
senting the preconditioned dual-time algorithm and its 
implementation in the diagonalized framework of the 
OVERFLOW code. Following this, we present several 
computational examples to verify the accuracy and ef- 
ficiency gains in various flow regimes of interest. The 
example problems include both simple test problems, 
such as the propagation of a Lamb vortex and Stokes’ 
first problem, as well as more practical situations, such 
as a round-jet in a low-Mach number cross-flow with 
ground-effect. The latter problem is representative of 
the Harrier flow fields mentioned earlier. Finally, we 
summarize the current status of the unsteady model- 
ing capability in OVERFLOW and point out the areas 
for future research and development. 


Preconditioned Dual-Time Algorithm 

Equations of Motion 

The Euler equations in generalized coordinates can 
be written as follows: 


dQ dE dF_ dG 

~di + di + d v + oq 


(i) 


Here, 
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and 


where J 1 = x^y n z^ -h x^y^z n + x^y^z^ x^y^z 7} 
x v i j n Z£ and U, V, and W are the contra- 
variant velocities. 

To carry out the iterations at each physical time 
step, an artificial time(r) term is explicitly introduced 


dQ_ dQ , dE dF dG 

dr dt <K dp dQ 


( 2 ) 


Convergence of the pseudo-time(sub-iterations) at 
each physical time step is important for obtaining an 
accurate transient solution. To optimize the perfor- 
mance of the pseudo-iterations, the pseudo-time term 
is written in terms of the primitive variable vector 
Q p = (P,u,v,v),T)/J and the preconditioning matrix 


r p- 

r dQp dQ dE OF dG _ 

r ”-dT + aT + ^ + W + sc 


(3) 


Here V p takes the form discussed in Ref. [9] and is 
given in the Appendix. 


Discretization 

Discretizing Eqn. 3 with first order finite difference 
for artificial time and second order backward difference 
for the physical time terms results in 

Q* +1 - Q p 3Q k + l -4Q n + Q n ~ l 

Ar + 2A t (4) 

+^£ fc+1 + S„F k+l + S c G k+1 = 0 
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Here, k is the pseudo-iteration counter, n is the time 
step counter and S represents spatial differences in the 
direction indicated by the subscript. After lineariza- 
tion, 


AQp 3(0* + r e Ag p ) - 4 Q" + Q’ 1 " 1 
At + 2At ^ 

+ S^E 11 + A k AQ p ) + 6,(F k + B k AQ p ) [0> 

+ 6 c (G k + C k AQ p ) = 0 


where A Q p = Q k+1 - Q k , A p , B p , and C p are the flux 
Jacobians with respect to Q p , and T e = Jqj- 

Rewriting this equation such that all terms evalu- 
ated at sub-iteration k or time steps n and n - 1 are 
on the right hand side and all terms multiplying AQ P 
are on the left hand side 


3 At 

Tp + Te 2Ai 


where 


+ A T(A$6t + B*6 n +Cfc) 


AQ P 


-- R k 

( 6 ) 


R k = -At 



4 Q n + Q n ~ l 
2 At 


+ S ( E k + 5 v F k + 6 c G k 


(7) 


The key step in the derivation of a diagonalized scheme 
rests in combining the pseudo- physical time-derivative 
terms on the left hand side into a single matrix. Ac- 
cordingly, we define S p = Ip + f^Te (see Appendix), 


(S p + ArAfe + A rB k 5 n + A rCfa) A Q p = R k 

(8) 

Multiplying through by r e and converting back to the 
conservative system, we get 

(l + ArT e SZ l A k 6 i + AtT'S! l B k S n 

V . (9) 

+ ATT e S; 1 C k S ( jAQ = r e S~ l R k 

where ,4, B , and C are the flux Jacobians with respect 
to Q . 

Approximate Factorization and 
Diagonalization 

Applying approximate factorization, 

(/ + ArTeS-M*^) (/ + Arr e 5 p - 1 H*<5 I) ) 

(/ + ArTeS; 1 ^^) A Q = T e S; l R k 

Now, r e S~ l A has the same eigenvalues as S~ 1 A p 
and the eigenvectors of T e S p x A are the columns of 
r e X£ where are the eigenvectors of S p l A p . Thus, 

r e S~ l A = T e S~ l A p T- 1 = (11) 


The scheme can be diagonalized using Eqn. 11 and 
takes the form 

T e Xz (^I 4- feAr^Ac j X' l X v (V + bArSjJV^ X v 1 A £ 

(i + 6Ar5 c A c ) X~ 1 T~ 1 AQ = 

(12) 

where £ = 1 + I ■ The eigenvalues A are discussed 
in the Appendix. 

Artificial Dissipation 

The diagonalized approximate factorization scheme 
employed in OVERFLOW uses central differences to 
discretize the spacial directions. It has been shown 
by many researchers that, to maintain uniform accu- 
racy across different flow regimes, it is important to 
formulate the artificial dissipation terms for the pre- 
conditioned system. In the dual-time context, this 
means that the artificial dissipation terms must be con- 
structed using the pseudo-time derivatives rather than 
the physical time- derivatives. Accordingly, a typical 
second order dissipation term can be written as 



where a(A) denotes the spectral radius of the matrix 

A. 

In the implicit scheme, the second-order dissipation 
terms appear in the implicit LHS operator as well. for 
the purposes of diagonalization, however, it is neces- 
sary to modify the spectral radius term a{T^ l A p ) with 
rr(5" 1 v4 p ). This approximation appears only on the 
LHS and has been shown by Buelow et al. [16] to have 
little or no effecct on the convergence efficiency. We 
further point out that the modification does not im- 
pact the accuracy of the discrete formulation. 

Preconditioner Definition 

The definition of the preconditioning matrix (r p ) 
is determined by the parameter M p (see Eqn. 17 in 
the Appendix). In the present case, this parameter is 
defined as 


M p 2 = M IN (M AX (Mf , Ml-Ml J, 1) (14) 

where M t is the local Mach number scale convention- 
ally used in steady preconditioning [17] and M 2 lin is 
a cutoff value usually set to 3 x M*,. M u is a new 
unsteady scale given by 


M u 


Lu 

7rA ta 


(15) 


where L u is a user specified unsteady length scale [16], 
and a is the speed of sound. Note that when M p = 1, 
preconditioning is not applied. 
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Note that the unsteady preconditioning varies the 
choice of the preconditioning parameter between the 
inviscid choice and the no preconditioning choice. For 
acoustic problems, the physical time-step size is small 
since it is based upon the acoustic CFL number. The 
unsteady preconditioning parameter then approaches 
unity (equivalent to the no preconditioning choice). 
On the other hand, for vortex propagation problems, 
the physical time step is large as it is selected based 
upon the particle CFL number. The corresponding un- 
steady preconditioning parameter then approaches the 
Mach number (equivalent to the steady precondition- 
ing choice). For intermediate physical time step sizes, 
the unsteady choice provides preconditioning values 
between these two limits. As discussed in Ref. [10], 
the unsteady preconditioning choice optimizes conver- 
gence efficiency of the pseudo-time iterations regard- 
less of the choice of physical time step. In most cases, 
it also insures accurate scaling of the artificial dissi- 
pation terms. The exceptions to this statement are 
discussed further in the Results section. 

Results 

Vortex propagation 

Vortex propagation is used as the first example of an 
unsteady flow field. Trends for this problem have been 
previously established in Ref. [10] based on stability 
analysis. Specifically, a Lamb vortex is defined with 
the following velocity distribution. 

V r (r,9) = 0 ; Vo = T ( 16 ) 

where F is the vortex strength and 0 is a character- 
istic radius. The initial conditions are specified using 
this velocity distribution and the vortex is propagated 
using the dual-time-stepping scheme described in the 
previous section. The present study is conducted at 
free stream Mach numbers of 0.1 and 0.001. To study 
the performance of the scheme with and without pre- 
conditioning, the sub-iteration convergence is studied 
for several time steps. The time step is characterized 
by CFL U = where u is the speed of propagation, 
At is the physical time step and Ax is the grid size. 
Ideally, for efficient time-marching of vortex propaga- 
tion problems, this CFL number should be of order 
unity. 

Here as well as in the other test cases presented 
in this paper, we compare the convergence of the 
sub-iterations at each physical time-step for three dif- 
ferent choices of the preconditioning matrix, namely, 
no preconditioning, steady preconditioning and un- 
steady preconditioning. The first choice corresponds 
to the traditional time-marching algorithm with dual 


time-stepping, which is the base scheme in the OVER- 
FLOW code. The second choice adopts the standard 
steady state preonditioning choice used in steady low’ 
Mach computations. Finally, the third choice concerns 
the present implementation and is given in Eqn. 14. 
We note that the length scale, L u , is typically selected 
to be the height of the domain, which is in accordance 
with the suggested value in Ref. [10]. 

Figure 1 shows the results for Af = 0.1 for four val- 
ues of CFL U . For the smaller values of CFL U (0.1 and 
1), the no preconditioning and unsteady precondition- 
ing cases converge rapidly, taking about 200 iterations 
to reach machine zero. On the other hand, the steady 
preconditioning choice is observed to yield poor per- 
formance. At the higher CFL numbers (10 and 100), 
the steady and unsteady preconditioning choices con- 
verge well, taking about 300 to 400 iterations to reach 
machine zero, w’hile the no preconditioning choice now r 
yields poorer performance. It is clear from these re- 
sults that the unsteady preconditioning choice essen- 
tially behaves like the no preconditioning choice for 
small time steps and like the steady preconditioning 
choice for large time steps. Importantly, the unsteady 
choice yields the best performance for all physical time 
step sizes. We point out however that for efficient 
and accurate solutions, the optimal physical time step 
choice would correspond to CFL U — 1. 

The advantage of using the unsteady preconditioner 
is observed to be more significant at even lowor speeds 
as Fig. 2 show's for a free stream Mach number of 
0.001. At this speed, for low values of CFL U the 
steady preconditioner does not converge in the sub- 
iteration. The unsteady preconditioner is observed to 
follow T the nonpreconditioned case as expected (see Fig. 
2(a)). At CFL U = 0.1, however, the unsteady pre- 
conditioner outperforms the non-preconditioned case. 
At CFL U = 1, the steady preconditioner is conver- 
gent, but like the nonpreconditioned case it displays 
poor performance. Again the unsteady preconditioner 
outperforms the other choices. At an even higher 
CFL U = 10 the steady preconditioner performs as well 
as its unsteady counterpart. It is thus clear that in all 
cases the unsteady preconditioning provides optimal 
performance of the pseudo-time convergence. Note 
that these trends arc similar to those presented in Ref. 
[ 10 ]. 

The subiteration convergence is important to study 
because it has an impact on the accuracy of the solu- 
tion. If the subiterations are not properly converged, 
the solution will be inaccurate. Of course, practical 
computations do not require that the sub-iterations be 
converged to machine zero. Therefore, to establish the 
amount of work required to obtain a solution, the num- 
ber of sub-iterations required for 4-order convergence 
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Fig. 1 Comparisons of residual convergence with and 

is plotted in Fig. 3 against CFL U . For smaller values 
of CFL U where the unsteady preconditioner performs 
as well as the no preconditioning choice, a relatively 
small number of sub-iterations are required. However, 
at these small CFL U , a large number of time steps 
is needed to march the solution a specified amount of 
time. Accordingly, the CPU time required for march- 
ing the vortex one rotation is plotted in Fig. 4. Here 
at low values of CFL U , the CPU time required is pro- 
hibitively high. For CFL U = 0(1), where the vortex 
is marching approximately 1 grid point per time step, 
the unsteady preconditioner requires the least number 
of sub-iterations and makes one rotation in the least 
amount of CPU time. At larger values of CFL U , the 
solution is no longer correct as the time step is simply 
too high to accurately capture the vortex convection. 
From these results, it is clear that the preconditioned 
dual-time scheme provides substantial savings in CPU 


without preconditioner for a Lamb vortex at A1 0.1 

time for low Mach flowfields. 

We next turn our attention to computational ac- 
curacy. Vortex computations are notoriously difficult 
and the second-order schemes considered here are espe- 
cially prone to excessive numerical damping. However, 
we are primarily concerned with a comparative assess- 
ment of the different preconditioning choices. Figure 5 
shows vorticity contours of the initial condition (Fig. 
5a) and the solution after one rotation through the (pe- 
riodic) domain for the various preconditioning choices 
(Figs. 5b, c and d). The non-preconditioned result 
(Fig. 5b) is completely wrong, with the vortex in the 
wrong location and its strength significantly damped 
out. Note that this is true even though the inner it- 
erations are fully converged. Thus, the discrepancy is 
due to the poor accuracy of the numerical discretiza- 
tion and not due to lack of convergence. In contrast, 
the steady (Fig. 5c) and unsteady preconditioning 
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M=0.001, CFL u =l 



Fig. 2 Comparisons of residual convergence with 

M = 0.001 

(Fig. 5d) results predict the correct location of the 
vortex and a more respectable vortex strength. These 
results clearly point to the benefit of preconditioning 
in limiting excessive numerical dissipation and hence 
preserving the accuracy of the numerical formulation 
as well as computational efficiency. 

Stokes’ first problem 

For the next test case, we consider the solution of an 
impulsively started plate, which involves the unsteady 
development of a viscous boundary layer. Commonly 
referred to as Stokes’ first problem, there is a closed 
form self-similar solution, which makes this an ideal 
test case for verifying both convergence and accuracy 
properties of the preconditioned scheme considered 
here. 

To systematically study the computational results. 


M=0.001, CKL U =10 



and without preconditioner for a Lamb vortex at 

we use an isotropic mesh with a fixed aspect ratio of 
unity. We use 11 grid points in the self-similar direc- 
tion and 101 grid points in the wall-normal direction. 
The domain height is 10 -4 for a A y = 10 6 . On this 
uniform Cartesian grid, the initial boundary layer is 
selected to extend to approximately 51 points off the 
surface. The initial condition is obtained from the 
analytical solution and the computations are carried 
out to follow the continued evolution of the boundary 
layer. 

We mentioned earlier that the characteristic length 
(L u ) that appears in the definition of the unsteady pre- 
conditioning parameter (Eqn. 15) was typically chosen 
to be the length or width of the computational do- 
main. Here, the relevant length scale is the height of 
the domain, which is given by 10 -4 . Figure 6 shows 
the sensitivity of the sub-iteration convergence to the 
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Lamb Vortex, M— 0,001 



Fig. 3 Number of sub-iterations required for 4 
order convergence for a Lamb vortex at A7 = 0.001 


Lamb Vortex, M=0.00! 



Fig* 4 CPU time required for 4 order 
converged (sub- iter at ion) Lamb vortex to make one 
rotation at M = 0.001 

selection of this parameter for a von Neumann number 
of 1.67 (VNN = = 1.67). It is clear that values 

of 10 -4 (and less) yield the best convergence rates. 
Figure 6 also shows the convergence rates with no pre- 
conditioning and the steady preconditioning choices* 
Indeed, it is clear that these choices yield substantially 
poorer sub- iteration convergence. 

To accurately capture the boundary layer growth, 
an appropriate physical time step for the numerical 
marching process must be chosen. The viscous time 
step or Von Neumann number (VNN) is the char- 
acteristic time-step of interest in this problem. It is 
desirable to use a von Neumann number of 1 which cor- 
responds approximately to a boundary layer growth of 
1 grid point per time step. To assess the convergence 
behavior of the low Mach preconditioner with respect 
to VNN, the number of sub-iterations required to con- 


verge the problem 4 orders of magnitude is plotted for 
a given time step in Fig. 7. For small values of V N N , 
the no preconditioning choice shows good convergence 
behavior. However, for higher values of VNN, the 
convergence is extremely slow and the solution can not 
be obtained in a reasonable number of sub-iterations. 
Using steady preconditioning, the convergence hangs 
before getting to 4 orders for VNN < 0.5. For higher 
values of VNN proper convergence can be obtained, 
but very large time steps can lead to poor resolution of 
the transient evolution of the solution. The unsteady 
preconditioner exhibits good sub-iteration convergence 
behavior regardless of the choice of V N A . 

The CPU time required to get four orders of sub- 
iteration convergence is shown in Fig. 8. It is clear 
that the unsteady preconditioning formulation yields 
the best CPU times for all choices of VNN. As men- 
tioned earlier, optimal efficiency and accuracy is ob- 
tained for VNN of order unity. Under these conditions, 
the no preconditioning case does not provide good con- 
vergence, while the steady and unsteady precondition- 
ing choices provide good sub-iteration convergence. 
Indeed, Figure 7 indicates that the dual time precon- 
ditioned formulation provides significant CPU savings 
over the traditional (noil-preconditioned) formulation. 

We have mentioned earlier that the choice of precon- 
ditioning also influences the accuracy of the computa- 
tions through the definition of the artificial dissipation 
terms. Table 1 compares the solutions obtained at a 
single point in the flowfield with the exact solution for 
several choices of VNN and preconditionings. It can 


Table 1 Solution comparison for Stokes’ first prob- 
lem with respect to various values of VNN 


VNN 

No 

Precond. 

Steady 

Precond. 

Unsteady 

Precond. 

Exact 

Solution 

0.017 

0.167 

1.67 

0.37151 

0.37177* 

0.37645* 

0.37617* 

0.37617* 

0.37618 

0.37465 

0.37606 

0.37619 

0.37622 

0.37622 

0.37622 


* Not converged 


be observed that the no preconditioning choice is not 
convergent for large VNN, while the steady precondi- 
tioning choice is not convergent for small VNN. The 
unsteady preconditioning, in contrast, is always con- 
vergent. The most accurate solutions are obtained by 
the steady and unsteady preconditioning choices when 
the VNN is order unity. It can also be observed that 
the no preconditioning result for small VNN is some- 
what inaccurate even though the solution is converged. 
This is because the poor scaling of the artificial dissi- 
pation terms yields a more diffusive formulation. For 
a more detailed explanantion, the reader is referred to 
Ref. [9]. Interestingly, the unsteady preconditioning 
result for small VNN is also somewhat in error. This 
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c) Steady Preconditioner 



Fig. 5 Comparison of solution with and without 


is because the unsteady preconditioning approaches 
the no preconditioning formulation at small VNN and 
the corresponding artificial dissipation terms become 
poorly scaled. This difficulty may be readily circum- 
vented by separating the preconditioning formulation 
used for the time-derivative (which is responsible for 
convergence effficiency) and that used for the artifi- 
cial dissipation terms (responsible for accuracy). Such 
a formulation may be implemented within a multiple 
pseudo-time framework as discussed in Ref. [18], but 
is beyond the scope of the present article. 

Cylindrical jet in cross flow with ground effect 

A complex unsteady flow field akin to the problem 
of a Harrier in ground effect is a jet in cross flow with 


preconditioner for a Lamb vortex at M — 0.001 

ground effect. This problem has been studied exper- 
imentally by Cimbala et. al. [19] with a jet at Mach 
0.13 issuing out of a 3 in. cylindrical tube in the cen- 
ter of a wind tunnel. The crosswind Mach number 
is 0.013. The jet subsequently impinges on the wind 
tunnel floor 9 in. away from the tube exit. 

For the present simulation a single grid with cylin- 
drical topology is used with stretching to assure 
proper viscous mesh spacing near the tube walls and 
ground. The mesh consists of approximately 250,000 
grid points. 

The experimental results of Cimbala ct al. [19] re- 
veal that the jet impinges on the ground and splays in 
all directions. The cross flow pushes the jet creating 
a horseshoe shaped vortex which surrounds the jet on 
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Comparison of subiteration performance with respect to changes in L u 
M =0.001 



Fig. 6 Parameter study to determine best L u for 
M = 0.001 on a grid of AR = 1 based on sub-iteration 
residual convergence of the Stokes’ first problem 


Stokes' first problem, M=0.001, AR=1 arid 



Fig. 7 Number of sub-iterations required 4 order 
convergence for Stokes’ first problem at M = 0.001 
on a grid of AR = 1 

3 sides. The ground vortex in turn is observed to pul- 
sate with a frequency of about 1-5 Hz, a result of the 
interaction between the jet with the low speed cross- 
flow. The general features of the flowfield are observed 
in Fig. 9, which shows a particle trace issued from the 
jet at a fixed point in time. The jet is seen to impinge 
the ground and a horseshoe vortex is seen in front of 
the jet. 

The computation of this flowfield with a traditional 
dual time scheme (without preconditioning) is stable 
only for small time steps. This stability restriction 
may be attributed to the poor convergence of the sub- 
iterations at large time step values. Our studies indi- 
cated that the maximum physical time step that could 
be reliably used was 5.5 x 10 -5 s. Using this time step), 


Stokes' first problem, M=0.001, AR=1 grid 



Fig. 8 CPU time required for 4 order 
converged (sub-iteration) Stokes’ first problem 
boundary layer to grow 10 grid points at M = 0 001 
on a grid of AR = 1 


more than 18,000 time steps are required to resolve a 
frequency of 1 Hz (that is characteristic of the ground 
vortex unsteadiness). Clearly, this would render the 
computation prohibitively expensive. Moreover, the 
non-preconditioned choice also leads to a more dissipa- 
tive numerical solution, which would also compromise 
the accuracy of the simulations. In contrast, when the 
preconditioned dual time scheme is employed, much 
larger time-steps may be stably employed arid reliable 
inner iteration convergence is obtained. The compu- 
tations presented in this paper were obtained with a 
physical time step size of 5.5 x 10 -3 s., which corre- 
sponds to about 180 time steps per 1 Hz cycle. We 
point out however that the larger physical time steps 
typically necessitate the use of more inner iterations 
(compared with five inner iterations used in the tradi- 
tional algorithm). In the present calculations, 50 inner 
iterations were employed at each physical time step. 
Thus, the preconditioned dual time scheme yields a 
potential CPU savings of about an order of magnitude. 
Figure 10 shows a plot of the L2 norm of the residual 
at each physical time step. The periodic unsteady na- 
ture of the flowfield is evident. The lowest frequency 
indicated in these results is approximately 1.8Hz. We 
point out that these results represent a preliminary 
effort of computing this complex flowfield. Detailed 
evaluations are presently underway and wall be the 
subject of a future article. Here, it is sufficient to state 
that the preconditioned dual-time formulation shows 
potential in computing complex multiple-timescale, 3- 
D problems with significantly improved accuracy and 
efficiency. 
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Fig. 9 The particle trace at an instance in ti 

Summary 

A preconditioned dual-time algorithm has been im- 
plemented in the OVERFLOW' code. The appioach 
follows the method of Buelow et al. to tailor the formu- 
lation within the diagonalized solution procedure used 
in the OVERFLOW' code. Modifications of the pre- 
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Fig. 10 The residual with respect to time step for 
the jet in cross flow 
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for a jet at A I — 0.13 in a Af — 0.013 cross flow 

conditioning selection to account for unsteady scales 
have also been included. The method significantly en- 
hances the capability' of the code, enabling it to handle 
unsteady flows over a wide range of Mach numbers and 
time scales. Specifically, the preconditioned dual-time 
method enhances accuracy by improving the scaling 
of the artifical dissipation terms used in the numer- 
ical discretization procedure and improves efficiency 
by optimizing the number of sub-iterations required 
at each physical time step. For low Mach number, low 
frequency problems, in particular, significant savings 
in CPU time are possible. 

The enhancements made possible by the method are 
demonstrated for both simple and practical flowfields. 
The simple test problems include the propagation of 
a Lamb vortex in an straight channel and the flow- 
field of an impulsively started plate. Computational 
results reveal that the dual-time scheme with unsteady 
preconditioning provides optimal sub-iteration conver- 
gence for all choices of time step size and Mach num- 
ber. Moreover, improved accuracy is also obtained 
with the method. 

Preliminary computations of a round jet in a low 
Mach crossflow* with ground effect have also been per- 

lautics and Astronautics 


formed. Here, the preconditioned dual-time scheme 
enables the use of much larger physical time step sizes, 
which is important for efficiently resolving the low fre- 
quency pulsations of the horshoe-shaped ground vor- 
tex. More detailed computations are currently being- 
performed and will be compared with experimental 
trends. The eventual goal is to apply the validated 
code to the computation of Harrier flowfields. 

The limitation of the present implementation is 
that the same preconditioning is employed for both 
the artificial dissipation formulation and for the time- 
derivative. The disadvantage of this is that the final 
solution is dependent upon the choice of time-step 
size. This difficulty may be effectively circumvented 
through using separate preconditioners for the dissi- 
pation and time-derivative terms. The former then 
controls the accuracy of the discrete formulation, while 
the latter controls convergence efficiency of the sub- 
iterations. Implementation of such a multi-level pre- 
conditioning formulation will be the subject of future 
work. 
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APPENDIX 

The Preconditioning Matrix 

The preconditioning matrix has the form defined 
in [16]. 


p'p 

0 

0 

0 

Pt 

upp 

p 

0 
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Upr 

vp'p 

0 

p 

0 

VpT 

wp’p 

0 

0 

p 

WpT 

(ph p + hop' p — 1) 

pu 

pv 

pw 

pflT + h,Qpr_ 


where p' p = — ^ and a is the speed of sound. Here e p = 

1Hy % M2 . However, in practice p[ p is never computed 
directly. Instead, computing the eigenvalues of S p l A p 
(Eqn. 22) involves the evaluation of the terms and 
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-jr. These terms are defined as follows. 
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Defining j3 l = 


Also, 
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(19) 


( 20 ) 


The values of 


A/p and 6 now controls the behavior of 
the preconditioner. The parameter b switches the be- 
havior of the preconditioner from unsteady to steady. 
For steady flows, 6=1 and /3' = M p . 

For unsteady flows, M p is defined by Eqn. 14. 

Definition of S v 
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( 22 ) 


Eigenvalues for S P l B p and S P 1 C P are similar. 

Eigenvectors of S p i A p 

For the repeated eigenvalues (A = bU), many choices 
of the eigenvectors are possible. The following choice 


keeps the eigenvectors associated with the linear eigen- 
values from becoming degenerate by assuring that the 
eigenvectors are non-vanishing independent of the ge- 
ometric orientation [17]. 
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The eigenvectors corresponding to the eigenvalues 
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For a perfect gas, the left and right eigenvector ma- 
trices can be written as follows. 
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and 


£ can be replaced with rj or £ to obtain the eigen- 
vector matrices for S^Bp, and S F l C p . 
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